C-------------------------------------------------------------------------------
C  y = Ax
C-------------------------------------------------------------------------------
      subroutine matvec_prod(Nx, Ny, Nz, A, x, y)
      implicit none
      
      integer Nx, Ny, Nz

      double precision A(Nx*Ny*Nz,4), x(Nx*Ny*Nz), y(Nx*Ny*Nz)

      integer NxNy, Nt, i

      NxNy = Nx*Ny
      Nt   = Nx*Ny*Nz

      do i=1,Nt
        y(i) = A(i,1)*x(i)
      end do
      
      do i=1,Nt-1
        y(i)   = y(i)   + A(i,2)*x(i+1)
        y(i+1) = y(i+1) + A(i,2)*x(i)
      end do
      
      do i=1,Nt-Nx
        y(i)    = y(i) +    A(i,3)*x(i+Nx)
        y(i+Nx) = y(i+Nx) + A(i,3)*x(i)
      end do
      
      do i=1,Nt-NxNy
        y(i)      = y(i)      + A(i,4)*x(i+NxNy)
        y(i+NxNy) = y(i+NxNy) + A(i,4)*x(i)
      end do
     
C      call output_cg(Nx,Ny,Nz, 'out.m', A, y, x)
C      stop
      
      return
      end
C-------------------------------------------------------------------------------
C  y = Ax
C-------------------------------------------------------------------------------
      subroutine matvec_prod_old(Nx, Ny, Nz, A, x, y)
      implicit none
      
      integer Nx, Ny, Nz

      double precision A(Nx*Ny*Nz,4), x(Nx*Ny*Nz), y(Nx*Ny*Nz)

      integer NxNy, Nt, i
      
      NxNy = Nx*Ny
      Nt   = Nx*Ny*Nz

      do i = 1, Nt
         y(i) = A(i,1)*x(i)
      end do

      i = 1
      y(i) = y(i)+A(i,2)*x(i+1)+A(i,3)*x(i+Nx)+A(i,4)*x(i+NxNy)

      do i = 2, Nx
         y(i) = y(i)+A(i,2)*x(i+1)+A(i,3)*x(i+Nx)+A(i,4)*x(i+NxNy)+
     +          A(i-1,2)*x(i-1)
      end do

      do i = Nx+1, NxNy
         y(i) = y(i)+A(i,2)*x(i+1)+A(i,3)*x(i+Nx)+A(i,4)*x(i+NxNy)+
     +          A(i-1,2)*x(i-1)+A(i-Nx,3)*x(i-Nx)
      end do

      do i = NxNy+1, Nt-NxNy
         y(i) = y(i)+A(i,2)*x(i+1)+A(i,3)*x(i+Nx)+A(i,4)*x(i+NxNy)+
     +          A(i-1,2)*x(i-1)+A(i-Nx,3)*x(i-Nx)+
     +          A(i-NxNy,4)*x(i-NxNy)
      end do

      do i = Nt-NxNy+1, Nt-Nx
         y(i) = y(i)+A(i,2)*x(i+1)+A(i,3)*x(i+Nx)+
     +          A(i-1,2)*x(i-1)+A(i-Nx,3)*x(i-Nx)+
     +          A(i-NxNy,4)*x(i-NxNy)
      end do

      do i = Nt-Nx+1, Nt-1
         y(i) = y(i)+A(i,2)*x(i+1)+
     +          A(i-1,2)*x(i-1)+A(i-Nx,3)*x(i-Nx)+
     +          A(i-NxNy,4)*x(i-NxNy)
      end do

      i = Nt
      y(i) = y(i)+A(i-1,2)*x(i-1)+A(i-Nx,3)*x(i-Nx)+
     +          A(i-NxNy,4)*x(i-NxNy)
     
      return
      end
C-------------------------------------------------------------------------------
C  z = (x,y)
C-------------------------------------------------------------------------------
      subroutine dot_prod(N, x, y, z)
      implicit none
      
      integer N
      
      double precision x(N), y(N), z
      
      integer i
      
      z = 0.
      do i=1,N  	      
        z = z + x(i) * y(i)    
      end do		      
      
      return
      end
C-------------------------------------------------------------------------------
C  y = x
C-------------------------------------------------------------------------------
      subroutine copy_vec(N, x, y)
      implicit none
      
      integer N
      
      double precision x(N), y(N)
      
      integer i
      
      do i=1,N  	      
        y(i) = x(i)
      end do		      
      
      return
      end
C-------------------------------------------------------------------------------
C  x(i) = alpha
C-------------------------------------------------------------------------------
      subroutine fill_vec(N, alpha, x)
      implicit none
      
      integer N
      
      double precision alpha, x(N)
      
      integer i
      
      do i=1,N  	      
        x(i) = alpha
      end do		      
      
      return
      end
C-------------------------------------------------------------------------------
C  z = x + alpha y
C-------------------------------------------------------------------------------
      subroutine lin_comb(N, x, y, alpha, z)
      implicit none
      
      integer N
      
      double precision x(N), y(N), z(N), alpha
      
      integer i
      
      do i=1,N  	      
        z(i) = x(i) + alpha * y(i)
      end do		      
      
      return
      end
